home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 039a / d3d.zip / BSPLSURF.C < prev    next >
Text File  |  1989-11-25  |  11KB  |  349 lines
  1. /* BSPLSURF: B-spline surface.
  2.    This program expects as D3D object file in which a grid of points is
  3.    defined.  The program asks the user to enter m and n: there must be 
  4.    exactly mn points in the file, numbered 1, 2,... mn.  It is assumed 
  5.    these points are arranged in a grid as follows:
  6.  
  7.          1        2        ...   m
  8.          m+1      m+2      ...   2m
  9.          2m+1     2m+2     ...   3m
  10.          .        .        ...   .
  11.          .        .        ...   .
  12.          .        .        ...   .
  13.          (n-1)m+1 (n-1)+2  ...   mn
  14.  
  15.    Neither m nor n must be less than 4.  The position of the points in
  16.    3-D space is free; no two points need have the same x-, y- or z- 
  17.    coordinates. If the object file contains a section beginning with 
  18.    "Faces:", then this section si ignored.  Two integera M and N are
  19.    to be entered; they denote numbers of intervals used in the surface-
  20.    fitting process.  There will be M intervals between points 1 and 2,
  21.    and N intervals between points 1 and m+1, and so on. The B-spline
  22.    curved surface that approximates the given points is written to an
  23.    output file in the D3D format.
  24. */
  25.  
  26. #include <stdio.h>
  27. #include <process.h>
  28. #include <alloc.h>
  29.  
  30. main()
  31. {
  32.    FILE *fp1, *fp2;
  33.    char str[30];
  34.    float *xx, *yy, *zz, x, y, z, v, u, u2, u3,
  35.       a00, a01, a02, a03,
  36.       a10, a11, a12, a13,
  37.       a20, a21, a22, a23,
  38.       a30, a31, a32, a33,
  39.       b00, b01, b02, b03,
  40.       b10, b11, b12, b13,
  41.       b20, b21, b22, b23,
  42.       b30, b31, b32, b33,
  43.       c00, c01, c02, c03,
  44.       c10, c11, c12, c13,
  45.       c20, c21, c22, c23,
  46.       c30, c31, c32, c33,
  47.       x00, x01, x02, x03,
  48.       x10, x11, x12, x13,
  49.       x20, x21, x22, x23,
  50.       x30, x31, x32, x33,
  51.       y00, y01, y02, y03,
  52.       y10, y11, y12, y13,
  53.       y20, y21, y22, y23,
  54.       y30, y31, y32, y33,
  55.       z00, z01, z02, z03,
  56.       z10, z11, z12, z13,
  57.       z20, z21, z22, z23,
  58.         z30, z31, z32, z33;
  59.  
  60.    int size, i, j, m, n, M, N, mn, l, k, I, J,
  61.       A, B, C, ii, jj, nn, mm, P,
  62.       h00, h01, h02, h03,
  63.       h10, h11, h12, h13,
  64.       h20, h21, h22, h23,
  65.         h30, h31, h32, h33;
  66.  
  67.    do {
  68.       printf("Enter m and n (both at least 4): ");
  69.       scanf("%d %d", &m, &n);
  70.    }
  71.    while (m < 4 || n < 4);
  72.  
  73.    printf("Input file:   "); 
  74.    scanf("%s", str);
  75.    if((fp1 = fopen(str, "r")) == NULL) { 
  76.       printf("No such file");
  77.       exit(1);
  78.    }
  79.    printf("Output file:  "); 
  80.    scanf("%s", str);
  81.    if((fp2 = fopen(str, "w")) == NULL) { 
  82.       printf("File Problem");
  83.       exit(1);
  84.    }
  85.  
  86.    mn = m*n;
  87.    size = (mn+1)*sizeof(float);
  88.    xx = malloc(size);
  89.    yy = malloc(size);
  90.    zz = malloc(size);
  91.    for(l=0; l<mn; l++) {
  92.       if(fscanf(fp1, "%d %f %f %f", &k, &x, &y, &z) != 4) {
  93.          printf("Wrong number of points in file");
  94.          exit(1);
  95.       }
  96.       if(k < 1 || k > mn) {
  97.          printf("Point number %d incorrect", k);
  98.          exit(1);
  99.       }
  100.       xx[k] = x;
  101.       yy[k] = y;
  102.       zz[k] = z;
  103.    }
  104.    fclose(fp1);
  105.  
  106.    printf("Enter M and N: ");
  107.    scanf("%d %d", &M, &N);
  108. /* Let us imagine we have n rows of m points, which gives us 
  109.    (n-1) * (m-1) rectangles.  Among these, only (n-3) * (m-3) 
  110.    rectangles will be approximated.  We consider each of these 
  111.    rectangles to be divided into M x N tiny rectangles, which 
  112.    we call 'elements'.  There will be A points in the output 
  113.    file for each horizontal grid line:
  114. */
  115.    A = (m-3) * M + 1;
  116.    printf("  i   j\n");
  117.    for(i=2; i<=n-2; i++)
  118.    for(j=2; j<=m-2; j++){
  119.       h11 = (i - 1) * m + j;
  120.       h10 = h11-1; h12 = h10+2; h13 = h10+3;
  121.       h00 = h10-m; h01 = h00+1; h02 = h00+2; h03 = h00+3;
  122.       h20 = h10+m; h21 = h20+1; h22 = h20+2; h23 = h20+3;
  123.       h30 = h20+m; h31 = h30+1; h32 = h30+2; h33 = h30+3;
  124.  
  125.       printf("%3d %3d\n",  i, j); /* Show we are busy */
  126.      
  127.       /* Here the 16 points are numbered as follows:
  128.  
  129.       P00  P01  P02  P03
  130.       P10  P11  P12  P13
  131.       P20  P21  P22  P23
  132.       P30  P31  P32  P33
  133.  
  134.       The inner region, within P11, P12, P21, P22, will be
  135.       approximated in this step (with the current values 
  136.       of i and j). */
  137.  
  138.       x00 = xx[h00]; y00 = yy[h00]; z00 = zz[h00];
  139.       x01 = xx[h01]; y01 = yy[h01]; z01 = zz[h01];
  140.       x02 = xx[h02]; y02 = yy[h02]; z02 = zz[h02];
  141.       x03 = xx[h03]; y03 = yy[h03]; z03 = zz[h03];
  142.  
  143.       x10 = xx[h10]; y10 = yy[h10]; z10 = zz[h10];
  144.       x11 = xx[h11]; y11 = yy[h11]; z11 = zz[h11];
  145.       x12 = xx[h12]; y12 = yy[h12]; z12 = zz[h12];
  146.       x13 = xx[h13]; y13 = yy[h13]; z13 = zz[h13];
  147.  
  148.       x20 = xx[h20]; y20 = yy[h20]; z20 = zz[h20];
  149.       x21 = xx[h21]; y21 = yy[h21]; z21 = zz[h21];
  150.       x22 = xx[h22]; y22 = yy[h22]; z22 = zz[h22];
  151.       x23 = xx[h23]; y23 = yy[h23]; z23 = zz[h23];
  152.  
  153.       x30 = xx[h30]; y30 = yy[h30]; z30 = zz[h30];
  154.       x31 = xx[h31]; y31 = yy[h31]; z31 = zz[h31];
  155.       x32 = xx[h32]; y32 = yy[h32]; z32 = zz[h32];
  156.       x33 = xx[h33]; y33 = yy[h33]; z33 = zz[h33];
  157.  
  158.       a33 = (x00-x03-x30+x33
  159.             + 3*(-x01+x02-x10+x13+x20-x23+x31-x32)
  160.             + 9*(x11-x12-x21+x22))/36;
  161.  
  162.       a32 = (-x00-x02+x30+x32 + 2*(x01-x31)
  163.             + 3*(x10+x12-x20-x22) + 6*(-x11+x21))/12;
  164.  
  165.       a31 = (x00-x02-x30+x32 + 3*(-x10+x12+x20-x22))/12;
  166.  
  167.       a30 = (-x00-x02+x30+x32 + 3*(x10+x12-x20-x22)
  168.             + 4*(-x01+x31) + 12*(x11-x21))/36;
  169.  
  170.       a23 = (-x00+x03-x20+x23 + 2*(x10-x13)
  171.             + 3*(x01-x02+x21-x22) + 6*(x12-x11))/12;
  172.       
  173.       a22 = (x00+x02+x20+x22 + 2*(-x01-x10-x12-x21))/4 + x11;
  174.  
  175.       a21 = (x02-x00-x20+x22 + 2*(x10-x12))/4;
  176.  
  177.       a20 = (x00+x02+x20+x22 - 2*(x10+x12)
  178.             + 4* (x01+x21) - 8*x11)/12;
  179.  
  180.       a13 = (x00-x03-x20+x23 + 3*(x02-x01+x21-x22))/12;
  181.  
  182.       a12 = (-x00-x02+x20+x22 + 2*(x01-x21))/4;
  183.  
  184.       a11 = (x00-x02-x20+x22)/4;
  185.  
  186.       a10 = (-x00-x02+x20+x22 + 4*(-x01+x21))/12;
  187.       
  188.         a03 = (x03-x00-x20+x23 + 3*(x01-x02+x21-x22)
  189.                 + 4*(x13-x10) + 12*(x11-x12))/36;
  190.  
  191.       a02 = (x00+x02+x20+x22 - 2*(x01+x21)
  192.             + 4*(x10+x12) - 8*x11)/12;
  193.  
  194.       a01 = (x02-x00-x20+x22 + 4*(x12-x10))/12;
  195.  
  196.       a00 = (x00+x02+x20+x22 + 4*(x01+x10+x12+x21)
  197.             + 16*x11)/36;
  198.  
  199.       b33 = (y00-y03-y30+y33
  200.             + 3*(-y01+y02-y10+y13+y20-y23+y31-y32)
  201.             + 9*(y11-y12-y21+y22))/36;
  202.  
  203.       b32 = (-y00-y02+y30+y32 + 2*(y01-y31)
  204.             + 3*(y10+y12-y20-y22) + 6*(-y11+y21))/12;
  205.  
  206.       b31 = (y00-y02-y30+y32 + 3*(-y10+y12+y20-y22))/12;
  207.  
  208.         b30 = (-y00-y02+y30+y32 + 3*(y10+y12-y20-y22)
  209.             + 4*(-y01+y31) + 12*(y11-y21))/36;
  210.  
  211.       b23 = (-y00+y03-y20+y23 + 2*(y10-y13)
  212.             + 3*(y01-y02+y21-y22) + 6*(y12-y11))/12;
  213.       
  214.       b22 = (y00+y02+y20+y22 + 2*(-y01-y10-y12-y21))/4 + y11;
  215.  
  216.       b21 = (y02-y00-y20+y22 + 2*(y10-y12))/4;
  217.  
  218.       b20 = (y00+y02+y20+y22 - 2*(y10+y12)
  219.             + 4* (y01+y21) - 8*y11)/12;
  220.  
  221.       b13 = (y00-y03-y20+y23 + 3*(y02-y01+y21-y22))/12;
  222.  
  223.       b12 = (-y00-y02+y20+y22 + 2*(y01-y21))/4;
  224.  
  225.       b11 = (y00-y02-y20+y22)/4;
  226.  
  227.       b10 = (-y00-y02+y20+y22 + 4*(-y01+y21))/12;
  228.       
  229.         b03 = (y03-y00-y20+y23 + 3*(y01-y02+y21-y22)
  230.                 + 4*(y13-y10) + 12*(y11-y12))/36;
  231.  
  232.       b02 = (y00+y02+y20+y22 - 2*(y01+y21)
  233.             + 4*(y10+y12) - 8*y11)/12;
  234.  
  235.       b01 = (y02-y00-y20+y22 + 4*(y12-y10))/12;
  236.  
  237.       b00 = (y00+y02+y20+y22 + 4*(y01+y10+y12+y21)
  238.             + 16*y11)/36;
  239.  
  240.       c33 = (z00-z03-z30+z33
  241.             + 3*(-z01+z02-z10+z13+z20-z23+z31-z32)
  242.             + 9*(z11-z12-z21+z22))/36;
  243.  
  244.       c32 = (-z00-z02+z30+z32 + 2*(z01-z31)
  245.             + 3*(z10+z12-z20-z22) + 6*(-z11+z21))/12;
  246.  
  247.       c31 = (z00-z02-z30+z32 + 3*(-z10+z12+z20-z22))/12;
  248.  
  249.         c30 = (-z00-z02+z30+z32 + 3*(z10+z12-z20-z22)
  250.             + 4*(-z01+z31) + 12*(z11-z21))/36;
  251.  
  252.       c23 = (-z00+z03-z20+z23 + 2*(z10-z13)
  253.             + 3*(z01-z02+z21-z22) + 6*(z12-z11))/12;
  254.  
  255.       c22 = (z00+z02+z20+z22 + 2*(-z01-z10-z12-z21))/4 + z11;
  256.  
  257.       c21 = (z02-z00-z20+z22 + 2*(z10-z12))/4;
  258.  
  259.       c20 = (z00+z02+z20+z22 - 2*(z10+z12)
  260.             + 4* (z01+z21) - 8*z11)/12;
  261.  
  262.       c13 = (z00-z03-z20+z23 + 3*(z02-z01+z21-z22))/12;
  263.  
  264.       c12 = (-z00-z02+z20+z22 + 2*(z01-z21))/4;
  265.  
  266.       c11 = (z00-z02-z20+z22)/4;
  267.  
  268.       c10 = (-z00-z02+z20+z22 + 4*(-z01+z21))/12;
  269.  
  270.         c03 = (z03-z00-z20+z23 + 3*(z01-z02+z21-z22)
  271.                 + 4*(z13-z10) + 12*(z11-z12))/36;
  272.  
  273.       c02 = (z00+z02+z20+z22 - 2*(z01+z21)
  274.             + 4*(z10+z12) - 8*z11)/12;
  275.  
  276.       c01 = (z02-z00-z20+z22 + 4*(z12-z10))/12;
  277.  
  278.       c00 = (z00+z02+z20+z22 + 4*(z01+z10+z12+z21)
  279.             + 16*z11)/36;
  280.  
  281.       /* In the output file we will be using point numbers k,
  282.          which are to be computed from the control variables 
  283.          i, j, I, J. In each rectangle, determined by their 
  284.          pair (i,j) and to be divided into N x M elements,
  285.          we use I, counting up to N, and J counting up to M.
  286.          To avoid duplicated points, I and J normally start
  287.          at 1, except for the lowest values (2) of i and j:
  288.          then I and J start at 0. (Of all rectangles that are
  289.          to be divided into elements, the one with i=j=2 is
  290.          the leftmost rectangle at the top.) If I=0 (which 
  291.          implies i=2) then k = (j-2_M + J + 1, otherwise
  292.          k = A + (i-2)NA + (I-1)A + (j-2)M + J + 1. This means
  293.          that in either case we have:
  294.  
  295.                k = {(i-2)N + I}A + {(j-2)M + 1} + J
  296.  
  297.          A is the total number of mesh points on a horizontal row.
  298.        */
  299.  
  300.        B = (j-2)*M + 1;
  301.        for(I=(i == 2 ? 0 : 1); I<=N; I++) {
  302.          u  = (float)I/N;
  303.          u2 = u*u;
  304.          u3 = u2*u;
  305.          C = ((i-2)*N+I)*A + B;
  306.          for(J=(j ==2 ? 0 : 1); J<=M; J++) {
  307.             k = C + J;
  308.             v = (float)J/M;
  309.             x =    ((a03*v + a02)*v + a01)*v + a00 +
  310.                      u*(((a13*v + a12)*v + a11)*v + a10)+
  311.                     u2*(((a23*v + a22)*v + a21)*v + a20)+
  312.                u3*(((a33*v + a32)*v + a31)*v + a30);
  313.  
  314.             y =    ((b03*v + b02)*v + b01)*v + b00 +
  315.                      u*(((b13*v + b12)*v + b11)*v + b10)+
  316.                     u2*(((b23*v + b22)*v + b21)*v + b20)+
  317.                u3*(((b33*v + b32)*v + b31)*v + b30);
  318.  
  319.             z =    ((c03*v + c02)*v + c01)*v + c00 +
  320.                      u*(((c13*v + c12)*v + c11)*v + c10)+
  321.                     u2*(((c23*v + c22)*v + c21)*v + c20)+
  322.                u3*(((c33*v + c32)*v + c31)*v + c30);          
  323.             fprintf(fp2, "%d %f %f %f\n", k, x, y, z);
  324.          }
  325.       }
  326.    }
  327.    fprintf(fp2, "Faces:\n");
  328.    nn = (n-3)*N;
  329.    mm = A-1;      /* mm = (m-3)*M */
  330.    for(ii=0; ii<nn; ii++)
  331.    for(jj=0; jj<mm; jj++){
  332.       P = ii*A + jj + 1;
  333.       /* We have to divide each (nearly rectangular) surface
  334.          element into two traingles to insure that the points
  335.          that will be used as vertices of polygons really lie
  336.          in the same plane.  The minus sign in, for example,
  337.          P -Q R will prevent D3D from drawing line segment PQ.
  338.          Since either face of each triangle will be visible,
  339.          we specify its vertices both counter-clockwise and 
  340.          clockwise:
  341.       */
  342.       fprintf(fp2, "%d %d %d.\n", P, -(P+A+1), P+1);
  343.       fprintf(fp2, "%d %d %d.\n", P+A+1, -P, P+A);
  344.       fprintf(fp2, "%d %d %d.\n", P, -(P+A+1), P+A);
  345.       fprintf(fp2, "%d %d %d.\n", P+A+1, -P, P+1);
  346.    }
  347.     fclose(fp2);
  348. }
  349.